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Abstract: A fifth-order relaxation scheme for the shallow water equations is presented in this paper. 
The scheme is based upon the fifth-order weighted essentially nonoscillatory (WENO) re- 
construction and the implicit-explicit Runge-Kutta scheme. The resulting scheme does not 
require Riemann solvers and the computation of Jacobians, so it has the advantages of 
relaxation schemes. The resulting method is applied to simulate the one-dimensional dam- 
break problems on flat and non-flat topography. The results demonstrate the robustness 
and effectiveness of the present method. The effect of bottom friction is also discussed. 
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1 Introduction 


The shallow water equations arise in many different applications, such as dam-break flows, 
prediction of water pollution, hydraulic jump, tidal flows in estuary, flood waves, etc. Recently, 
high-resolution relaxation schemes originally suggested in [1] have attracted much attention. 
In comparison with upwind schemes such as the Godunov scheme, relaxation schemes do not 
require the it Riemann solvers and the computation of Jacobians. These features make the 
relaxation schemes particularly suitable for those systems where the Riemann problem is dif- 
ficult to solve or when it is not possible to perform analytical expression for Jacobians. The 
relaxation schemes were successfully implemented to the incompressible Euler equations and 
Hamilton-Jacobi equations. 

As for application and implementation in solving the shallow water equations, Delis and 
Katsaounis!) had adopted the first-order and second-order relaxation schemes introduced in 
[2]. Recently, Seaid!3) presented a general framework to generalize the relaxation schemes to 
higher orders of accuracy and developed a third-order scheme for the shallow water equations. 
A fourth-order relaxation scheme for hyperbolic conservation laws can be found in [4]. In this 
paper we present a fifth-order relaxation scheme for approximating solutions of one-dimensional 
shallow water equations. Our method is a high-order extension of the scheme in [1,3,4]. In 
application, the present scheme is used to simulate dam-break problems. 
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2 Governing equation 


The one-dimensional shallow water equations can be written in conservation form as 


U: + F(U); = S(U) , (1) 
where 
h h 0 
U= | F(U)= E , S(U)= 
hu hu? + igh? —ghB, — ghSfex 


Here, h is the water depth, u is the velocity, B(x) is the bottom topography, g is the gravitational 
constant and Sfr = n2ululh—3 denotes the friction slope, where n is the roughness coefficient. 


3 A fifth-order relaxation scheme 


Based on the ideas presented in [1,2], the one-dimensional shallow water equations (1) can 


be replaced by the following relaxation system 


BU: OV oa 

Ot Oz (2) 
Ov ðU 1 

a a a 


where V € R?, 7 > 0 is the relaxation rate and A = diag{a;, a2} is a positive diagonal matrix. 
Under certain conditions the solution of (2) converges to the solution of the original problem 
(1) as r — 0. Without loss of generality, let us consider the following uniform spatial grid 
. An Al 
zj =jAT, Tj = G Ë 5) Ae. 


3 


Introduce cell average of the variable U in I; = |z;_ poja ] as 


U;(t) = > I U(a, t) dz. 


The approximate point value of U at (x;,1,¢) is denoted by U;42 (t). A spatial discretization 
to (2) in conservation form is given by 


Ms 4 i Vay Vins)“ 800 : 
3 

OV; 1 1 

a + Ags ies —U;_1) = -7V3 — F(U);). 


Assuming that the cell-averages {U;} are known, our goal is to compute the cell point values 
U;+, and V;,1, the flux F(U); and the source S(U);. The computation proceeds in two 
steps. As in [3], the first step is to perform a reconstruction from the given cell-averages 


U(a) = 2 R;(2; U)x;(z), (4) 
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where x; is the characteristic function of the interval J; and R;(z;U) is a polynomial defined 
in I}. Given such a reconstruction, the point-values of U at the interface points are denoted by 


Uj. = R(x j+4;U) and Uj_, = R,(#;_3;U). 


To improve resolution and reduce numerical dissipation, the fifth-order WENO reconstruction!®! 
is utilized in this paper. With this background the computation of U;, ı from fifth-order WENO 
2 


reconstruction is 


where U, ith and V; j+} are the k-th components of U; 44 and V} +45 respectively. The qí ı and 
2 


gE r) 


j, are 
Gea > SU; F 20541 = EUjsa, ay = -5 j-1 + ŠU; + Ujan, 
aa EU; = TU + EU; qs = Tu; = Uj $ “Ujta, 

g, = Uj $ U; = Use, a”, = -Uj k SU; + SU; 
The weights w, and ®r are given by 
. . S; _ d ; 
uea a= a a Se ae a= Ay i =0,1,2 


Here dọ = da = 3, dı = d = 3, dg = dọ = b and € = 107°. The smoothness indicators, ĝi, 


are calculated by 


13 1 

Bo = pU ~ 2Uj41 +Uj+2) + racy — 4U;41 + Uj+2)°, 
13 1 

A= iz U- — 2U; + Uj41)” + zU —Uj41)’; 


13 2 1 2 
Bo = jz Ui-2 —2Uj;-1+ U;) + zU- — 4U;-1 + 3U;) j 
Then, the k-th components of variables V + A2U is discretized by 
(V+ vaU); = Rj(£j41; V + VaxU), (V- VarU) 5,3 = Rys1 (25433 V — vaU), 
where V is the k-th components of V. The point values are obtained by 


Via = NANE (25435 V + vVakU) — Ry41 (25435 V — VaxU)), 


= z (Ri(£j} p; V + VaeU) + Rii (2j+3; V — VaxU)). 
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Here Uj} and V;,1 are the k-th components of U j+4 and V,,1, respectively. The next step 
is to approximate the flux F(U); and the source s(U); j by the Sipen’ s quadrature. 
To implement the time discretization, the formulation (3) is rewritten in common ordinary 


differential equations form 


dY 
Z =LY)- Y), (5) 
where 
U; S(U); — D: V. i 0 
Y= , L(Y)= (Usea , L(Y) = 
V; —AD,V; V; — F(U) 
Here 1 
DV: = Zz (Vir — Vi-4). 


The third-order implicit-explicit Runge-Kutta scheme proposed in [6] is utilized in this paper. 
When applied to system (5) it can be written as 


l-1 
Ki =Y" + At Ý GimL(Km = D CimL(Km), l= 1,2,3, 


m=1 


yrtl y” + ats dL(Ki) -Ê a= alta. 
t=1 


The 3 x 3 matrices C = (Gim) and C = (cim) are given by 


0 0 0 0 0 0 
C= y 0 0j, C=)10 y D4; 
y-1 2-2y 0 0 2~2y y 


where y = 2+y3 . The coefficient vectors d = (d,,dz,d3)7 and d = (d;,d2,d3)" are given as 
= 1 1\T 1 1\T 
d=(0,5,5 3 d=(0,5,5) . 


4 Numerical examples 


In this section, several test problems are presented. The results demonstrate the performance 
of the proposed scheme. In all the numerical examples below, the computational parameters 
are 7 = 1075, a = sup|u — vgh| and /az = sup|u + vghl. 

Example 1: One-dimensional dam-break problem 

This problem is used to test the accuracy of the algorithm. Consider a horizontal rectangular 
channel with 2000m long. A dam is located at 1000m. The initial upstream and downstream 
water depths, h, and ha, are given by following two cases: 1) hy = 10m,hg = 5m; 2) 
hy = 10m, ha = 0.5m. Thus the water depth ratios ha/hu are 0.5 and 0.05. The computational 
domain is discretized into 200 cells. The bed slope is set to zero. Result at t = 50s for depth 
ratio 0.5 is shown in Figure 1. It can be seen that the shock and rarefaction waves are well 
resolved by the present method. Figure 2 shows the result at t = 50s for depth ratio 0.05. In 
this case, the roughness coefficient is considered and set as 0, 0.02 and 0.04, respectively. 
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Figure 1: Water height, ha/hu = 0.5 Figure 2: Water height, ha/hu = 0.05 


Example 2: Dam-break problem over a continuous riverbed 

We verify our scheme on the dam-break problem over a continuous riverbed. This test case 
is taken from [7]. The bottom function, B(x) = 1.398 — 0.347 tanh(8x — 4). The computaitonal 
domain is (0,1]. The initial conditions are 


1.0 x< 0.6, 
h(x,0) = hu(z,0) = 0.0. 
0.2 z¢>0.6, 


The gravitational constant is set to g = 1.0. A uniform grid with 100 grid points is chose in 
this case, and the numerical solution is plotted together with the reference solution which is 
computed by the present scheme with 3000 grid points. Figure 3 shows the result at t = 0.25. 
The simulated result demonstrates the effectiveness of our method. 


Figure 3: Water height and discharge plot for dam-break problem over a continuous riverbed 
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